
#	Andy Harris
#	Date: 20 Aug 2014
#	Purpose: Function to implement MC simulations.

#	libraries
library(foreach)
library(doMC)
registerDoMC(12)
#	Functions
xsumx <- function(x){x/sum(x)}
sampfunc <- function(rnx, x, tmpNN, tmpparams){sample(rnx, size = round(tmpNN*tmpparams), replace = T, prob = x)}
source("nameEst.R")
source("nameEstW.R")
source("mc.table.R")

#	Start Code

start.time <- Sys.time()
#Parameters from which to simulate names data: group proportions we attempt to estimate.
param.list <- list(
c(rep(1/6, 6)),
c(0.30, 0.30, 0, 0, 0.1, 0.30),
c(0.55, 0.25, 0.05, 0.05, 0, 0.1),
c(0.05, 0.8, 0.05, 0.025, 0.05, 0.025),
c(0.8, 0.05, 0.05, 0, 0.05, 0.05),
c(0.01, 0.11, 0.02, 0, 0.01, 0.85)
)


#Sample sizes: number of names in each simulation.
nn.list <- c(2000, 4000, 8000, 16000, 32000)

#Number of simulations
nsims <- 5000
cen <- readRDS("freqsUS.Rdata")

train <- foreach(ii = 1:ncol(cen), .combine = cbind) %dopar% {
	(cen[,ii]/sum(cen[,ii]))
}
rownames(train) <- row.names(cen)
colnames(train) <- names(cen)

#################
#United States
#################

tab.ne.us <- mc.table(nsims = nsims, param.list = param.list, nn.list = nn.list, traindata = train, namefunc = nameEst, classify = NULL)

tab.new.us <- mc.table(nsims = nsims, param.list = param.list, nn.list = nn.list, traindata = train, namefunc = nameEstW, classify = NULL)

classifyUnif <- train/rowSums(train)
classifyUS <- cen/rowSums(cen)
tab.unif.us <- mc.table(nsims = nsims, param.list = param.list, nn.list = nn.list, traindata = train, namefunc = NULL, classify = classifyUnif)

tab.pop.us <- mc.table(nsims = nsims, param.list = param.list, nn.list = nn.list, traindata = train, namefunc = NULL, classify = classifyUS)





#################
#Kenya
#################

train <- readRDS("nameConditionals.Rdata")
param.list <- list(
	rep(1/ncol(train), ncol(train)),
	c(0.4, 0.1, 0.2, 0.05, 0.05, 0.05, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025),
	c(0.05, 0.1, 0.4, 0.05, 0.05, 0.20, 0.025, 0.025, 0.025, 0.025, 0.025, 0.025),
	c(0.1, 0.25, 0.25, 0.05, 0.05, 0.25, 0.025, 0.025, 0,0,0,0),
	c(0.1, 0, 0.1, 0.05, 0.25, 0.15, 0, 0.1, 0, 0, 0.25, 0),
	c(0.125, 0.025, 0.1, 0.05, 0.3, 0.3, 0, 0.1, 0, 0, 0, 0),
	c(0.025, 0.025, 0.8, 0.05, 0, 0, 0, 0.1, 0, 0, 0, 0),
	c(0.8, 0.025, 0.025, 0.05, 0, 0, 0, 0.1, 0, 0, 0, 0)
)
nn.list <- c(2000, 4000, 8000, 16000, 32000)

tab.ne.ke <- mc.table(nsims = nsims, param.list = param.list, nn.list = nn.list, traindata = train, namefunc = nameEst, classify = NULL)
tab.new.ke <- mc.table(nsims = nsims, param.list = param.list, nn.list = nn.list, traindata = train, namefunc = nameEstW, classify = NULL)

classifyUnif <- train/rowSums(train)
tab.ke <- mc.table(nsims = nsims, param.list = param.list, nn.list = nn.list, traindata = train, namefunc = NULL, classify = classifyUnif)




#Make main tables in article.

#table 1, panel 1:
t1p1 <- tab.ne.us$point/tab.unif.us$point

#table 1, panel 2:
t1p2 <- tab.ne.us$point/tab.pop.us$point

#table 1, panel 3:
t1p3 <- tab.new.us$point/tab.ne.us$point

#table 1, panel 4:
t1p4 <- tab.ne.ke$point/tab.ke$point

#These are the four panels in table 1 of the article. Assembled by hand into latex table.
save(t1p1, t1p2, t1p3, t1p4, file = 'mctablemat.Rdata')


end.time <- Sys.time()

start.time
end.time

